"""
list oriented functions for computation
of statistical information.
"""

import math

def mean(list, start=None, stop=None, bad=-999.0):
    length = len(list)
    if not start: start = 0
    if not stop: stop = length
    sum = wgt = 0
    for i in range(start, stop):
	datum = list[i]
	if datum == bad: continue
	sum = sum + datum
	wgt = wgt + 1
    if wgt == 0: return bad
    return sum / wgt

def anom(list, mean, bad=-999.0):
    length = len(list)
    if mean == bad: return [bad] * length
    result = [bad] * length
    for i in range(length):
	datum = list[i]
	if datum == bad: continue
	result[i] = datum - mean
    return result

def sigma(list, mean, start=None, stop=None, bad=-999.0):
    length = len(list)
    if mean == bad: return bad
    if not start: start = 0
    if not stop: stop = length
    s2 = n = 0
    for i in range(start, stop):
	datum = list[i]
	if datum == bad: continue
	n = n + 1
	s2 = s2 + pow(datum - mean, 2)
    if n == 0: return bad
    variance = s2 / n
    sigma = pow(variance, 0.5)
    return sigma

def meanAndSigma(list, start=None, stop=None, bad=-999.0):
    length = len(list)
    if not start: start = 0
    if not stop: stop = length
    sos = sum = wgt = 0
    for i in range(start, stop):
	x = list[i]
        if x == bad: continue
        sos = sos + x * x
        sum = sum + x
        wgt = wgt + 1
    if wgt == 0:
        return bad, bad
    else:
        mean = sum / wgt
        sigma = math.sqrt(sos / wgt - mean * mean)  
        return mean, sigma

def normalize(list, mean, sigma, bad=-999.0):
    if mean == bad: return [bad] * len(list)
    def f(x, mean=mean, sigma=sigma, bad=bad):
	if x == bad: return bad
	return (x - mean) / sigma
    return map(f, list)

def anomaly(list, mean, bad=-999.0):
    if mean == bad: return [bad] * len(list)
    def f(x, mean=mean, bad=bad):
        if x == bad: return bad
        else: return x - mean
    return map(f, list)

def trend(list, bad=-999.0):
    length = len(list)
    n = 0
    Sx = Sx2 = Sxy = Sy = 0.0
    for i in range(length):
        y = list[i]
	x = float(i)
        if y == bad: continue
        n = n + 1
        Sx = Sx + x
        Sx2 = Sx2 + x*x
        Sxy = Sxy + x*y
        Sy = Sy + y
    if n == 0: return (bad, bad)
    m = Sxy - (Sx / n) * Sy
    m = m / (Sx2 - (Sx / n) * Sx)
    b = (Sy - Sx * m) / n
    return (m, b)

def change(list, trend=None, bad=-999.0, offset=1):
    if trend == None:
	trend = trend(list, bad)
    m, b = trend
    return m * (len(list) - offset)
